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ABSTRACT 

We consider semi-discrete Galerkin approximation schemes in connection 
with inverse problems for the estimation of spatially varying coefficients and 
boundary condition parameters in second order hyperbolic systems typical of 
those arising in 1-D surface seismic problems. Spline based algorithms are 
proposed for which theoretical convergence results along with a representative 
sample of numerical findings are given. 
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1. Introduction. In this paper we consider computational techniques for 
the following class of inverse problems: For the system 

2 

(1.1) p(x)^ = -^(E(x)^) t>0, 0<x<l , 

9t 3X 3X 

(1.2) k^v(t,0) = s(t;k) 

(1.3) ||(t,l) + k2|J(t,l) = 0 

(1.4) v(0,x) = <f>(x) , v^(0,x) = ij;(x) , 

given observations {y,-,-} for {v(t. ,x.)}, choose. from some admissible set. 

"best" estimates for the parameters p, E, k^, k 2 » k. These problems are 
motivated by certain versions of the so-called "1-D Seismic Inversion Problem" 
(see, e.g. [1 ], [8]). Roughly speaking, one has an elastic medium (e.g., 
the earth) with density p and elastic modulus E. A perturbation of the system 
(explosions, or vibrating loads from specially designed trucks) near 
the surface (x=0) produces a source s for particle disturbances v that travel 
as elastic waves, being partially reflected due to the inhomogeneous nature 
of the medium. An important but difficult problem involves using the observed 
disturbances at the surface or at points along a "bore hole" to determine 
properties (represented by parameters in the system) of the medium. In the 
highly idealized 1-D "surface seismic" problem, one assumes that data are 
collected at the same point (x=0) where the original disturbance or "source" 
is located. In addition to this hypothesis which cannot be true, other unreal- 
istic special assumptions are made about the nature of the traveling and re- 
flected waves. Although the standard 1-D formulations are far from reality. 
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exploration seismologists have developed techniques for processing actual field 
data (performing a series of experiments and "stacking" the data) so that the 
1-D problems are generally accepted as useful and worthy subjects of investi- 
gation. Consequently, numerous papers (for some interesting references, see 
the bibliographies of [1 ], [8]) on the 1-D problems can be found in the research 
literature. 

In many formulations of the seismic inverse problem, the medium is assumed 
to be the half-line x>0 (with x=0 the surface) while in others (especially 
some of those dealingwith computational schemes) one finds the assumption of an 
artificial finite boundary (say at x=l) at which no downgoing waves are reflected 
(an "absorbing" boundary). While there are several ways to approximate such a 
condition in 2 or 3 dimensional problems (see [12], [21]), for the 1-D formulation 
this condition is embodied in a simple boundary condition of the form (1.3); 
here V .2 ~ ^ E(l)/p(l) and one can view this boundary condition as resulting 
from factoring the wave equation (1.1) at x= l and imposing the condition of 
"no upgoing waves" at x=l. 

Equation (1.1) is a 1-D version of the equations for an isotropic elastic 
medium while (1.2) represents an "elastic" boundary condition at the surface 
x=0 (k-| represents an elastic modulus for the restoring force produced by the 
medium) . 

As is the case in many inverse or "identification" problems, the problems 
described above tend to be ill-posed (including a computationally undesirable 
instability) unless careful restrictions are imposed on the admissible parameter 
class (for some discussions of these aspects, see [1 ], [10]). We shall not 
focus on this aspect here. Rather, the purpose of our presentation in this 
paper is to demonstrate the feasibility of a certain theoretical approach and 
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certain approximations in developing computational schemes for problems in 
which there are i) unknown boundary parameters and ii) unknown spatially 
varying coefficients in the system equations. 

We choose the "1-D seismic inverse problem" involving (1.1) - (1.4) as 
a test example to exhibit the efficacy of our ideas. However the technical 
features and notions we present are of importance in a number of other ap- 
plications. There are rather easily motivated and fundamental problems in 
dealing with large elastic structures (large space structures - e.g. beam- 
like structures with tip bodies) that involve estimation of boundary condition 
parameters. In these cases the models are often hybrid models with distributed 
system (Euler-Bernoulli , Timoshenko) state equations and ordinary differential 
equation boundary conditions (see, for example, [2], [9], [18], [20]). A 
second class of problems for which the techniques introduced in this paper 
have immediate use are related to bioturbation [7], [13], This is the mixing 
of lake and deep-sea sediments by burrowing activities of organisms. Under- 
standing of this phenomenon is fundamental to geologists in interpreting geologic 
records contained in sediment core samples. The best models to date involve 
parabolic state equations (for a nonuniform "mixing chamber") with unknown 
parameters in the boundary conditions describing the flux into and out of the 
chamber. 

In our approach here we employ the Trotter-Kato theorem to obtain theoretical 
convergence results (assuming regularity of parameter sets to guarantee 
existence of solutions to the inverse problems) for spline approximation 
schemes for the states. Boundary parameter estimation is treated directly 
via mappings that iteratively change the parameter-dependent spline basis 
elements into "conforming" elements (i.e., elements which satisfy the appropriate 
boundary conditions). We deal only with estimation of regular spatially- 
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varying coefficients in (.1.1), where again splines are used for parameters 
in a secondary approximation. Estimation of discontinuous coefficients 
(including location of the discontinuities) in problems such as those that 
are the focus of our attention in this paper can be effectively treated 
theoretically and numerically in a framework similar to that here using, 
for example, tau-Legendre state approximation schemes [4 ]. 

We turn then to the estimation problem for (1.1)- (1.4). It is theoretically 
and numerically advantageous to deal with homogeneous boundary conditions by 
transforming the problem so that the source term s in (1.2) appears in the 
initial data and in a term in the state equation. We make the transformation 
u = V + G where (here represents differentiation with respect to t) 

G(t,x;q) = - (^)s(t;k) + ( ^ )x2(x-l )s(t;k) 
and obtain the system 

u^(t.O) + q3u(t,0) = 0 

(1.5) 

u^(t,l) + q^u^(t,l) = 0 

u(0,x) = i(x;q) , u^(0,x) = i(x;q). 

Here the forcing function F is given by 

F(t,x;q) = q^(x) {-(i)s(t;k) + (^)x2(x-l )T (t;k)} 

H 3 q3q4 

3 4 
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where here and throughout we adopt the notation q = (q-j* ^ 2 * ^3’ ^41 

q^ = p, qg = E. ‘’s '^r ^4 ^ 4* transformed initial conditions have 

the form 

Hxiq) = <^(x) - (^)s(0;k) + (^)x^(x-l )s(0:k) 
i(x;q) • m - (^)s(0;k) + (^)x^{x-l)V(0;k) . 

We assume henceforth that we have observations = ^^11* ***’ ^im^ * 

i = l, 2, ..., n, corresponding to w(tj*,q) = (u(t^-,x-|), ...» u(t. ,X|^)) where 

* ^ 

u is the solution of (1.5). For a criterion in determining a best estimate q 
of the parameters we use a least-squares function 

(1.6) J(q) = Z ly,- - w(t.;q)l^ 

i=l ^ ^ 

which we seek to minimize as q ranges over some admissible parameter set Q. 

We remark that in the event our observations n^* = *'*’ ’^im 

original system (1.1)-(1.4), we may apply directly the theory and techniques 

of this paper by considering in place of (1.6) the criterion 

(1.7) J(q) = Z In. + G(t*;q) - w(t.;q)l 

i=l ^ ^ 

where 6(t^iq) = (G(t^jXijq), ...« 6(t« tX^^^jq)) . 

We make some standing assumptions to facilitate consideration of our 

problem in subsequent discussions. We shall search for q in a set 

Q c= C(0,1) X H^(0,1) X R X R X (we shall sometimes write Q as x x 

1 2+k 

Q, X Q- X Q_). We further assume that Q is compact in the C x H x R 

o ^ 0 

topology, and that there exist positive constants 
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q, . q, 

, i=l, 2, 3, 4 such 

that 



^i 1 1 q,- 

for q . e Q . , 

1 = 1, 2, 


^3 1 -'»3 i ^3 

for q 3 e Q3 , 

and 


q4 

foi" q4 e Q 4 . 


Finally, 

we assume ((le H^(0,1) , 

ij>eH^(0,l) , and 

s(* ;k) e H^(0,T) for each 


keQs* where e [0,T], T<®, and that k -> s(*;k) is a continuous mapping from 
Qg to H^(0,T). 

We turn next to the theoretical foundations of the approximation schemes 
we propose to use in solving our inverse problem of minimizing J over Q, 
subject to (1.5). 
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2 , Abstract Formulation. 

The object in this section is to lay the theoretical foundation for the 
problem. First, we shall write our partial differential equation as an abstract 
ordinary differential equation in a Hilbert space, then determine a set of 
approximating ordinary differential equations. Each of these abstract equations 
will have an associated identification problem; the original will be referred 
to as (ID), the approximating problem will be referred to as (ID^). We 
shall use the theory of semigroups to obtain existence and uniqueness of 
solutions to the differential equations. We can then fit our problem into the 
theoretical framework developed in [ 5 ], and deduce that, under conditions 
stated there (reiterated below for clarity), one can solve (ID ) for each N, 
and these parameter estimates thus obtained will "lead to" a solution of (ID). 

The equation (1.5) can be rewritten as a first order system, motivating 
the use of a product (V(q) x L (q)) of two spaces to be our Hilbert space X(q). 

Define V(q) to be H^(0,1) with inner product defined by = 

P 2 DvDwdx - q2(0)q2v(0)w(0). (D denotes the spatial differentiation operator 
— ). It can be readily shown that for any qeQ, V(q) is a Hilbert space, and 
moreover, the assumptions made about Q imply that the V(q) norm is uniformly 

equivalent to the H^ norm as q ranges over Q. Let Vg(q) contain those elements 

of V(q) which satisfy the elastic boundary condition* i.e., Vg(q) = 
{veV(q)nH^(0,l)|Dv(0) + q3v(0) = 0}. 

We define L (q) to be H (0,1) with inner product given by <v,w>q ^ = 

1 2 
q^vwdx , and note that for each qeQ, L (q) is a Hilbert space and its norm 

is uniformly equivalent to the standard H^ norm as q ranges over Q. 

2 

As described earlier, we take X(q) = V(q) x L (q) with inner product 

given by <x,y>^ = <x-i»yT>\/(q) + ^^2’^2^0,q ^ “ (x^,X 2 )^ and y = (y^,y 2 )^). 

It is clear from our remarks above that for qeQ, X(q) is a Hilbert space, and the 
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X norm is uniformly equivalent to the norm as q ranges over Q. We can 

formally write (1.5) as an abstract equation in X(q): 


( 2 . 1 ) 


z(t) = A(q)z(t) + G(t;q) 


z(0) = Zg(q) 


/u(t,*) 

where we have identified z(t)eX(q) with ( 1 . The boundary conditions 

\u^(t,.); 

are incorporated into the domain of A(q) by defining domA(q) = {(y)eVg(q) x 
H^(0,l)|v(l) + q^Du(l) = 0}, and A is the unbounded linear operator given by 


A(q). = 


(l/q^)D(q2D) 0 


The function G and the initial condition are given by 


G(t;q) 




It can be shown that for each qeQ, A(q) is the infinitesimal generator of 
a CQ-semigroup, T(t;q) on X(q), so that we have the existence of mild solutions 
to (2.1), given by 

t 

(2.2) z(t;q) = T(t;q)z«(q) + / T(t-s;q)G(s;q)ds 

0 

with z(*;q)eC(0,T;X(q)). In this context, the inverse problem can be stated as 

(ID) Given observations y = {y^.}”_-|, minimize J(z(*;q).y) over qeQ 
subject to z(*;q) satisfying (2.2). 
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Here, J(q) a J(z(-;q),y) = “ 5(t^,q)|^ where 5(t^.q) = 

{z,(t.,x^;q), .... z^(t^.,x^;q)) and z^ denotes the first component of z. 

To prove that for each q, A(q) generates a CQ-semigroup, one can use the 
Lumer-Phillips Theorem ([15], p.l6). To employ this theorem, one must show the 
operator is dissipative, densely defined, and satisfies a certain range statement. 
To demonstrate the .dissipativity of A(q), we take fedomA(q), qeQ, and compute 
(with an integration by parts) 


<A(q)f,f>q 


{ " ) 

y\(l/qi)D(q2Dfi)/ 



= "■^2’^l"V(q) <n/qi)D(q2Dfi).f2>o.q 



= - q2(0)q3fi(0)f2(0) - q2(0)Df^ (0)f2(0) + q2(l )Dfi (1 )f2(l ) 


= - q2(l)q4(Dfi(l))^lO . 

By relating domA(q) to other subsets (see [Ifl for details) which are known to 
be dense in one can easily argue that for each qeQ, domA(q) is dense 

in X(q). One can also argue that R(x-A(q)) = X(q) for some x>0, by demon- 
strating that given/^^^e X(q), there exists ^^jedomA(q) such that 


Xu - V 

•(l/q-|)D(q2Du) + 

This is equivalent to solving the following two point boundary value problem: 
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- (Vq])D(q2Du) + X^u = Xf^ + 

Du(0) + q2u(0) = 0 
XU(1) + q^Dud) = f^O) 

for ueH^CO.l), and setting v(x) = xu(x) - f^(x). 

If we let y = u - (l/q 4 )x^(x-l)f^(l) the above problem is transformed 
to an equivalent one with homogeneous boundary conditions: 

(-l/q^)D(q2Dy) + x^y = F 

Dy(0) + q 3 y.( 0 ) = 0 

q 4 Dy(i) + Xy(l) = 0 
2 

where FeL (q). One can then use the theory of self-adjoint operators (again 

seeCl4) to argue that a solution exists for any FeL^(q). 

We now turn to the approximation of our equation (2.1). We shall obtain 
N 

a solution z to an approximating equation (to be discussed in detail below) 
in a finite dimensional subspace of X(q), denoted X*^(q). Specifically, let 
S (a ) represent the standard subspace of cubic splines corresponding to 
the partition " VN (see pp. 78-81 of [16]): then, 

given qeQ, we take X^(q) to be that subspace of S^(a^) x S^(a*^) whose elements 
satisfy the boundary conditions corresponding to q (i.e., X*^*(q) c domA(q) ). 

Let Bj , j = -1,...,N+1,. be the B-spline basis elements for S^(a*^). Then 
X^(q) is the (2N + 3)-dimensional subspace spanned by the following set of basis 
f uncti ons : 
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Let P^(q):X(q)-»-X^(q) denote the orthogonal projection of X(q) onto 
X^(q), i.e., given feX(q), p'^(q)f is that element in X^(q) which satisfies 
|P^(q)^'- fl 1 |g- all 9 a Pof" each qeQ, we define an operator 


I 
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A^(q) on ^(q) given by A^(q) = P^(q)A(q)p\q), and then the approximating 
equation to (2.1) is written as: 

z^(t) = A^(q)z'^(t) + P^(q)G(t;q) 

(2.3) 

2^(0) = P^(q)zQ(q) 

where z^(t)eX^(q). Using the fact that A(q) is closed, p'^(q) is bounded, 
and the Closed Graph Theorem, one finds that A^(q) is bounded. The operator 
A (q) inherits the dissipativity of A(q), and therefore it follows that for 
each qeQ, A^(q) is the infinitesimal generator of a CQ-semigroup of contractions 
T^(t;q) on X(q). It is readily seen that Tf^(t;q) leaves X^(q) invariant. Thus, 
for each qeQ and each N=l,2, ..., there exists a unique mild solution 
z*'*(* ;q) e C(0,T;X*^(q)) of (2.3), which can be expressed as 

(2.4) z^(t;q) = T^(t;q)P^^(q)z (q.) + /V(t-s;q)P^^(q)G(s;q)ds. 

0 

The associated approximate identification problem is given by 

(ID^) Given observations y = , minimize J(z^(-;q),j) over qeQ 

subject to z^(-;q) satisfying (2.4). 

Here, J^(q) = J(z^(-;q),y) = ^ ,|y. - ?^(t.,q)|^ where C^(t.,q) = 

i=l ' 

(z-j (t^. ,x^ ;q) , ..., Zi (t. ,Xjj^;q)) and z!j^ denotes the first component of z^. 

Since X^(q) is finite dimensional, (2.3) is in fact a system of 2N+3 
ordinary differential equations, which can be solved using standard numerical 
packages. Similarly, there are numerical packages available to solve (ID^), 
provided solutions exist and we have some computationally feasible representation 
for q^ and qg. A detailed description of our numerical implementation, 
including a discussion of possible representations of q^ and will be deferred 
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to subsequent sections. First, our concern is to determine under what conditions 

solutions of (ID^) exist and how they relate to a solution of (ID). This is 

the subject of the next theorem, a slight modification of that given in [5 , p. 820]. 

1 2+k 

Theorem 2.1. Assume Q is compact in the C x H x R topology. If 
q-»-ZQ(q), q-»-P^(q)f, q -»-T^(t;q)f , fe X = X(q) are continuous in this same 
Q-topology, with the latter uniformly in te[0,T], then 

(i) There exists for each N a solution q^ of (ID^) and the 

''N L 

sequence {q } possesses a convergent subsequence q ^q. 

(ii) If we further assume that, for any sequence {q"^} in Q with 
« • • 

q^-^q, we have Iz^(t;q^) - z(t;q)l --^0 as j uniformly in 

qj 

A 

te[0,T], then q is a solution of (ID). 

The reader may, at first glance, find the convergence statement of (ii) 
suspect in that (t;q'^ ) e X^'(q^) and z(t;q)eX(q) , but this statement is 
meaningful in view of the following observation. In defining the spaces 
V(q), L^(q), and X(q), it was noted that V(q), L^(q), and X(q) are uniformly 
equivalent to h\ H®, and H^xH®, respectively, as q ranges over Q. This 
implies that the X(q) are setwise equal as q ranges over Q. To be technically 
precise, we should use the canonical isomorphism when relating an element of 
X(q'^) to its counterpart in X(q), but to simplify our presentation, we shall 
throughout abuse notation and omit the isomorphism. 

It is easily seen from the form of ZQ(q) that q -*■ ZQ(q) is continuous. 

It is also true that for our p'^(q), T^(t;q) we have q P^(q)f and q -»■ T^(t;q)f 
continuous; this will be readily seen from the matrix representations for our 
approximating scheme, and so further discussion is postponed until Section 5. 
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The next theorem gives sufficient conditions for the hypothesis. of (ii) 
from Theorem 2.1 to hold. 

Theorem 2.2. Let q be arbitrary in Q such that q q as N ->• » (recall 

convergence is in the C x h’ x topology). Suppose that the projections 

P%) are such that l(P^(q^ -l)fl ^^0 as N -> » for all feX(q), that 

feX(q) implies |T^(t;q*'*)f - T(t;q)f | as uniformly in te[0,T], and 

that lzg(q ) - Zg(q)| 0 as N ^ Then the mild solutions z^(t;q*^) of 

q 

(2.3) converge to the mild solution z(t;q) of (2.1) uniformly in te[0,T]. 

The proof of this theorem, which is based on a standard "variation-of- 
constants" representation for solutions z and z^ in terms of the semigroups 
T and T , essentially follows immediately from Theorem 3.1 of [ 5 , p. 823]. 

One only needs to verify that our spaces, operators, etc. satisfy the conditions 
required in [ 5 ]. 

It is clear from the continuity of q ^ Zg(q) that |zg(q^) - Zg(q)| ^0 as q^^q 

qN 

It remains only to show the convergence of the projections and the semigroups. 

The main result of the next section is the convergence of the semigroups; the 
convergence of the projections is obtained as an intermediate proposition. 

In summary then, at the end of the next section, we will be able to deduce 
from Theorem 2.2 that z*'^(t;q*'^) converges to z(t;q) whenever q*^ -»■ q, and hence 
by Theorem 2.1 we are assured that the sequence of iterates {q*'*} we obtain 
by solving (ID^), has a subsequence which converges to a solution, q, of (ID). 
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3. Convergence Arguments. 

This section will be devoted to establishing the result; For each convergent 
sequence q in Q, and for any feX(.q), lT^Ct;q^)f ~ T(t;q)fl ^ 0 as 

N ^ uniformly in te[0,T3. As explained in the previous section, this 
convergence result is crucial in arguing that z^(t;q^) -»■ z(t;q) whenever 
q^ q, which in turn is necessary to ensure that our candidate (the limit of 
our approximating subsequence) is indeed a solution to our inverse problem. 

We shall first prove a slightly different form of convergence of the semi- 
groups using the following version of the Trotter-Kato Theorem [3 3, 

Theorem 3.1 . Let (B,l*l) and (B^,l*l^). N = 1, 2, ..., be Banach spaces and 
let : B->B^ be bounded linear operators. Further assume that T(t) and T^(t) 
are CQ-semigroups on B and b” with infinitesimal generators A and A , respectively. 
If 

(i) lim jn^'^fL = |f| for all feB, 

N->co 

(ii) there exist constants M, u independent of N such that 
lT^(t)lf^ 1 Me“^, for t > 0, 

(iii) there exists a set VCB, VO dom(A), with (Xq-A)I? = B for 
some Xq> 0, such that for all feV we have 

lA^f - n^Af luj -^0 as N ^ , 

then lT^(t)n^f - n^T(t)flf^ ^ 0 as N for all feB, uniformly in t on compact 
intervals in [0,«). 

N 

It will be a standing assumption throughout this section that q q 
in Q with this convergence in the C x x topology. Let B = X(q) with 
norm denoted by |-1^, B^ = X(q'^) with norm j • ] for N = 1 , 2, . . . , A = A(q) 
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with corresponding semigroup T(t) = T(t;q), and = A^(q^) = P^(q^)A(q^)P^(q^) 
with corresponding semigroup T^U) = T^Ct;q^) (as described in Section 2). For 

each N, : X(q) X(q^) will be a bounded linear operator which will map elements 
of domA(q) into elements of domA(q^). Define 


g^(x) = axpCCq^- q^jx) - (x^/ 2)[q2- q^JexpCq,- q-]; 


given f let be defined by 


n'f = 


9 ^1 

,(q4/q4)g% 


The functions g^ are defined so that as N -> «, g'^(x) 1, and D'^(g^(x)) *>0 

for any positive integer j, where in each case the convergence is uniform 
in x€ [0,1]. 

A simple computation demonstrates that if fedom A(q), then n^fedom A(q^). 
For each N, n is a bounded linear operator from X(q) to X(q*^), but moreover, 
the set of operators {n^} is uniformly bounded. This statement can be proved 
using the assumptions on Q and the properties of g*'* mentioned above. Similar 
comments apply to the proof of our first proposition. 


Proposition 3.1 . 


For any feX(q), |n'^f 


f I 0 as N 
cN 


In order to argue the convergence of the infinitesimal generators, we 
shall need error estimates for the spline approximations and their derivatives. 
These will be variations of estimates such as those found in [19], modified 
to take into account our q-dependent norm, and the presence of the operator 


17 


The following notation will be used throughout this section. Given a 
vector function f, we shall use f^ or (f)^ to denote the i component of f. 
Given the scalar function h, will denbte the standard cubic spline 
interpolant of h (thus l\ e S^Ca^)). For a vector function f =(fM » 
will be the vector whose components are the spline interpolants of the 

f I^f \ 

components of f, i.e., I*^f = i j and I^f e S^(a^) x s^(a^). The 

interpolant of f which satisfies the boundary conditions corresponding to 
q will be written as Ig(q)f. While I^f interpolates f^ and at the values 
^i/I^^_Q and the derivatives of f-| and f£ at 0 and 1, will interpolate 

f^ and fg at the values additionally satisfy 

CD(Ig(q)f)^.J(0) + q^(Ig(q)f)^J(0) = 0, or equivalently, 

CD(lg(q)f),O(0) = - qgf^.CO) for i = 1. 2, 
and 

[Clg(q)f) 2 ] 0 ) + q4[0(lB(‘l)f)l3O) = 0» o'" equivalently, 

[D(I^(q)f)l](D = - (1/q4)f2(l). 

We note that if f satisfies the boundary conditions involving q, then 

Ig{q)f = 

The first estimates involve cubic interpolants for scalar functions. 

Lemma 3.1 . If h e H^, then 

iD^(h-I^h) Iq ^ 0 as N -»■ » , 

10(h-I^h)lQ < N"^lD^(h-I^h)|Q<. N'^jD^hlg . 

Ih - i\1q < N'2lD2(h-I^h)lQ < N'^lD^hlg • 


The convergence statement of this lemma follows immediately from the 
density of in H^, the estimates of Theorem 6.9 of [19], and the first 
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integral relation (4.15) of [19j. The estimates follow from (4.24) and (4.25), 
respectively, of [19] and the first integral relation. 

One can use the results of Lemma 3.1 and the equivalence oftlie X(q) 
and x norms to derive similar statements for the interpolants in the 
X(q) norm". 

Lemma 3.2 . If f e x and q e QC C x x then 

|D{A-f)|^ < K^dO^Cfi - l"f,)|^+ |D2(f2- 
where K^, Kg are constants which are independent of f, q, and N. 

Again, due to the equivalence of norms, the Schmidt inequality of [13, Thm. 1.5] 
can be modified and used component-wise to give a Schmidt type inequality in 
the X(q) norm. 

Lemma 3.3. If f e S (a ) x S (a^) and qeQ, then |Df|q ^ KgNjfl , where 
Kg is a constant independent of f, N, and q. 


The preceding estimates can be used to establish convergence properties 
for the canonical projections P^(q'^) where q^ q in Q. 

Proposition 3.2. If feX(q), then 
lP^(q^)f- f 1 0 as N -> «• . 

q 

Proof . First consider f e domA(q)f1(H^ x H^). For such f, n'^’f e domA(q^)fl(H^ x H^) 
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and Ig(q*^)n*^f = We use Lemma 3.2 in the triangle inequalities below 

to derive 

|P^q^)f- f < |P^(q^)[f - n^fllqN + |P^q^)n^f - + In^f- fiqN 

< 2|n^f-flqN + llB(q^)n^f - 

= 2ln^f-f|qM + llVf - n^flqN 

< 2|n^f-f|qfl + K^N"''(lD^(n^f)^l^ + lD2(n^f)2lo)^^^ . 

Thus we have lP^(q^)f- f|qN bounded by terms which we can show converge to 

N 

zero using Proposition 3.1 and the properties of g . 

The P*^(q^) are uniformly bounded, and the set domA(q)n(H^ x H^) is dense in 
X(q), hence one can use standard arguments to conclude that the statement of 
the proposition holds for all feX(q). 

Proposition 3.3 . For each feX(q), I(P*^(q^) - On^flq^ ^ 0 as N and 
for each f e domA(q)fl(H^ x H^), lD[(P*'*(q*'^) - I)n^f]|q|^j ->■ 0 as N -> «. 

Proof . The first statement is proved within the proof of Proposition 3.2; 
specifically, it was shown that |P^(q^)n*^f - n^f|qf^ 1 \ lD^(n*'^f)l Ig 

+ |D^(n^f)2lg)^''^. 

The proof of the second statement is obtained from the following triangle 
inequality (here we also use Lemmas 3.3, 3.2); 

lD(P^(qN)nWn|qf^< lD(P^q^)n^f - l^(n^f)) + iDd'^(n^f) - n^f)|^f^ 

1 K3N|P^(q^)n^f - I^(n^f)|qN + I.D(i^(n^f) - n^f)lqifj 

1 K3N|(P^(q^)-I)nNf|qN + K3Nln^f - iVflq^ + |D[I^(n^V) - n^f]|qN 
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< 2K3N|lVf-n''f|^N * |DDVf-n”f]|^„ 

< (2K,K3-MC2)C|D^[(ll"f), - l\n"f),]|2+ lO^CCAjj - 

Thus the conclusion lD[(P**(q“)- I)n’*f]|^N * 0 as N -► » follows from the 
observation that for i = 1,2 

f |D^[l"fl - fj]lo + iD^Cf,- - (Alfllo 

< 2|D^C(n''f ), - fp|„ + |D^[l"f, - f,]|„ 

with the latter terms approaching zero because of the properties of 
and Lemma 3.1, respectively. 

In later arguments, it will be helpful to have bounds (in the and tP 

norms) on one component of an element of X in terms of a bound (in the X(q) norm) 

on the entire element. Thus, we consider for feX(q), |f|^ = If-jly^qy + 

|f 2 lo,q which is equivalent to iDf^l^ + lf-,|Q + lf 2 lQ , so that there exist 

constants and k 2 such that iDf^lp < k^|f|q and [fglg < k 2 |flq . Similarly, 

|Df|q = |D^ilv(q) + l°^2lo,q equivalent to |D^f^l^+ |Df^l^+ lDf 2 |Q 

so we infer the existence of constants k 2 and k^ such that iD^f-jlQ ± k^lDfjq 
2 2 

and lDf 2 lQ <. k^lDfjq. For future reference, we combine and label these 
observations as 

IWlloihlflq 

(3.1) |D^f,|^lk3lDf|2 

1 ksiflq * • 

It is now possible to state and prove the following convergence theorem. 
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Theorem 3.2 , Suppose q in Q (convergence is in the topology). 

Then 

lT^(t;q^)n^f - n^T(t;q)flqN -^0 as N , 
for all feX(q), uniformly in t on compact intervals in [0,»)* 


Proof . The result is an immediate consequence of Theorem 3.1, once the 
hypotheses of that theorem have been shown to hold. Part (i) follows from 
Proposition 3.1, while part (ii) holds since T^(t;q) and T(t;q) are 
contraction semigroups for each N and qeQ. It remains only to verify (iii), 

0 O 

for which we take V to be the set domA(q)n(H x H ). Let f eV. Then 

|A^(q^^)n^f-n^(q)flqH = |P^(q^)A(q^)P^q^)n^f - n^A(q)f 

1 lP^q^)[A(q'^)P^{q^)^^f-^^(q)f]lqN + lP^q^)n^A(q)f - n^(q)f I^n 

1 |A(q^)P^(q^)n^f-n^A(q)flqN + l(P^(q^) - I)n^A(q)f 


= e-,(N) + GgCN). 


It follows directly from Proposition 3.3 that e 2 (N) 0 as N -»■ « . We must 
work harder to establish that ei(N) ->■ 0. We begin by breaking the norm into 
its two components and treat each separately. Thus 


[e-|(N)r = 


Jl/q!!')D(qjD) 


(p*^(q^)n^f) 


(l/qi)D(q2D) 



- I (P (q )n f )2 “ g '^2'V(q‘^) 


+ I(l/q5')[>[q^0(p"(q'')n''f),] - 


Cq3/q4)g''O/q,)D[2i20f,]|2_^„ 
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2 [S,(N)]^ + [«2CN)]^ 
We first observe that 


«,(N) i |(pN(q'‘)nNf)2 - CqJ/S4)A2lv(qN) ^ ICCq^) ‘ ’]9%lv(qN) 
= l(p''(q'')ll"f)2 - (n"f)2l„(<^) + l((qS|/q4) - Dg^lvf,/!) • 


It is more convenient, and due to the equivalence of the norms, it is sufficient, 

to establish the convergence in the norm. This can easily be done for the 

first term by invoking Proposition 3.3 and the inequalities (3.1). An argument 

can be made for the second term based on the properties of the g*^ and the 
N 

convergence q -> q. 

We turn now to the estimation of Using the equivalence of the 

L^(q''*) and norms, and the inequalities (3.1), we establish the following 
chain of inequalities: 

N 

« 2 <N) = 1-^ OCq” 0(P''{q'')A),]- 3 ") ^ q'< 

3 i 34 qi 


which is equivalent to 




3NJ. 




nN 


N J_ 


Dq.Df 


2"'ll0 



°^2 ^4 N 
^1 ^4 


1 '0 
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1 D2(P^q")II''f), -^DVf),lo+l^ “Slo* 

qi ^1 ^1 ‘’4 


nnN 

Dq2 

' oN 

'll 


D(P^q^)n^f), D(n^f)T 1 q+ l^D(n^f)i g^)Df^ |, 

^1 ^1 ^4 


.N 


■»'}' -N 


< A:4lJD[(p''(q"H)n’'f]l,N* l^l„|D((p"!q")-i)A),L. 

9l ‘’l 


•v -N 


!| D2(n"f),-^(^)g" D^f 1„+ 1^ D(n"f), g")Of,lo 

h qi ^4 ‘’i ‘’4 


We thus see that SgCN) can be bounded by four terms which go to zero as N-»-«; 

the convergence of the first two terms is the result of Proposition 3.3 and 

the convergence of q^ to q, while the convergence of the second two can be 

N . N 

argued using the properties of g and q -»-q. 


N 

We can use this theorem, the convergence properties of the operators n 
(Proposition 3.1), and the semigroup properties of and T, to establish the 
final result we need, as a corollary. 

Corollary 3.1 . Suppose q*^-»-q. Then 

lT^''Ct;q^)f- T(t;q)flqN ->• 0 as N « 

for all feXCq), uniformly in t on compact intervals in [0,~)- 


24 


We can now invoke the results (see Theorems 2.1 and 2.2) stated in Section 2 
to conclude that q (.obtained there as the lijnit of an approximating subsequence, 
{q is a solution to the identification problem. 
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4. Parameter Approximation. 

N 

In Section 2,. we pose the problem of minimiziing J‘ (q) over Q. The 
arguments underlying Theorem 2.1 yield that (under certain assumptions) 
each (approximate) problem has a solution q^, and for any convergent 
subsequence {q^^}> with q^^ q» we have q is a solution of the original 
identification problem. Recall, however, that q^i and are functional 
coefficients, and hence each of the approximate optimization problems is in 
fact infinite dimensional in nature. In this section, we discuss some 
methods for approximating these infinite dimensional optimization problems 
by finite dimensional ones, thus providing numerically tractable problems. 

This, of course, results in a second, or parameter, approximation that must 
be considered. 

In Section 5, we shall present the results of several numerical test 
examples. To reduce ill-posedness (see the comments in Section 1) we set 
q^ = p E 1 and search for q 2 = q 3 » q 4 » i<» with q 2 the only functional 

unknown. We therefore restrict our theoretical discussions here to this 
case. (We note however that in principle, our methods and ideas can be applied 
to the estimation of both p and E.) 

An approach that one might take would be to assume ana priori parameter- 
ization for q 2 * Thus the estimation of the unknown function becomes the esti- 
mation of a set of unknown constants appearing in the parameterization. The 
convergence theory developed thus far is directly applicable to this method. 
However, it would only yield results for best approximates (through the cri- 
terion on state observations) to q 2 within the fixed a priori parameterization 
class . Little can be said about convergence to a "best fit parameter" q 2 from 
the original parameter set Q. 


r 
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An alternate approach, which does not require. qualitative (e.gf, shape) 
assumptions about the parameter class, is to search for the unknown parameter 
in a sequence of sets which are finite dimensional approximations to the 
set Q. For example, one might search for the unknown parameter in sequences 
of classes of linear combinations of spline (or members of any other suitably 
chosen approximation family) basis elements. 

We shall consider here two cases: as a set of linear spline inter- 

polants, and as a set of cubic spline interpolants. For both cases we need 
to generalize the theory developed in Section 2, since we now have a "double 
index" (reflecting approximations for both the parameter and the state space) 
sequence of iterates, which we would like to argue converges to a solution of 
the original identification problem. 

To be specific, let Q = Q2 Q3 x Q4 x Qg = x and assume we 

have a mapping i*^ : Qg h\ For I the identity map, define = i^ x (l)^'*'*^, 
i.e., for qeQ, we have l”(q) = (i^(q2), q^, q^). 

Let = T^(Q). We assume 

(4.1) The set (Q^)g h i*'^(Q2) is compact in H^. 

(4.2) For c|2 6Q2» i (^2^ ^ ^2 as M -> « , and this convergence is 

uniform in q2eQ2* 

The original set Q is assumed to 

be compact in x so it follows from (4.1), the definition of and 

Theorem 2.1 that for each N and M, a solution qj| exists to the problem of 

minimizing J*'* over Q^. From the definition = ^^(Q), we see that there exists 
-N 

q^ e Q such that I (q^) = qj^ for each N and M. But the compactness of the 
original set Q then implies the existence of some subsequence {qjjj} and an 
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_Ni 

element q eQ such that qn^'^ q in Q; moreover, this subsequence may be 

Ic <<v 

chosen so that both Nj » and Mj^ ». The limit q is in fact a solution to 

the problem of minimizing J over Q; this claim is verified as follows: From 

-Ni 

the definition qj^'' we have 




This implies 


(4.3) 1 0^^(I^'^(q)) 


for q e Q . 


."Ni - , _Ni, ,_Ni - , . -Ni - 

But - q I < |I (q^J) - q^^'^ [ + - q | , and thus q,^'’ q in Q 


Nj M|^ ->■ ® follows from (4.2), the definition of and qJJ^ ^ q . If 
take the limit in (4.3) as Nj, Mj^ -»• <», we see that J(q ) £ J(q) for qeQ. 
Here we have used Theorem 2.2 with the observation that the convergence 
statement z^(t;q*^) -► z(t;2f) for any q*^ ^ q is still valid if replaced by 
^*^(t;q'^) z(t;q) as j, N -»■ », for any q"^ ->■ q; this can be seen using a re- 
indexing argument. These remarks are summarized in the following theorem. 


-Ni- 


as 


we 


Theorem 4.1 . Let = I*'**(Q) where (4.1) and (4.2) are satisfied. Let qjjj| be 

N M 

a solution to the problem of minimizing J over Q . Then for any convergent 
-Ni -Ni 

subsequence (qj^ ) with Nj, ~ and ->■ q , the limit q is a solution to 
k k 

the problem of minimizing J over Q. 


We first consider the above results applied to the' case where the Q*'*' are 
sets of linear spline interpolants. Let (4^) represent the subspace of 
piecewise linear splines corresponding to the partition x^. = 
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and let ->■ S^(A^) denote the standard linear spline interpolating 

operator. If, in addition to assuming Qg is compact in we assume Q2 
satisfies Qg C. | |D^q2lQ 1 K} , then it is not difficult to show 

that (4.1) and (4.2) are true for and i^ as defined above. From a standard 
representation result for linear interpolating splines [19, p.l2], we infer the 
continuity of the operator i**^ as a mapping from to h\ and the compactness 
of (Q^)2 = i^(Q2) io follows immediately. To establish (4.2) we appeal to 
standard estimates such as (2.17) and (2.18) in [19]. Having verified (4.1) 
and (4.2), we now state 

1 2+k 

Theorem 4.2. Suppose Q = Q2 ^ Q3 ^ Q4 ^ Q5 is a compact subset of H x R 
with Q2 additionally satisfying Q2 C {q2 e j lD^q2lo — “ I^(Q) 

where = i^ x (I)^^*^, and i^.is the linear spline interpolating operator. 

If qjll represents a solution obtained from minimizing over then for 

•'Ni ''N ~ 

any subsequence of {qj^} such that as N^, M|^ q^.^ -> q in Q, we have 

^ k k 

that q is a minimi zer for J over Q. 

Under slightly stronger assumptions on the set Q, we can develop a similar 

3 M 

convergence result using cubic spline approximations to q2* Let S'’(a”) be 

the subspace of cubic splines corresponding to the partition A^, and let 

i^ : S^(a^) denote the standard cubic spline interpolating operator 

(see Sections 2 and 3 for details). We assume Q2 is a compact subset of 

0 0 

satisfying also Q2 «= {q2 6H |D q2lQ ± K}. We again may use standard 

M 

interpolating spline representations (see [19, p. 45] ) to conclude that i is a 
continuous operator from to H^, from whence it follows that (Q^)2 is 
compact in H^. To verify (4.2), we again refer to (4.19) and (4,20) in 
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[19]. Thus we have 

1 2 +k 

Theorem 4.3. Suppose Q = Qg ^ Q 3 Q 4 ^ Q 5 is a compact subset of C x R 
with Qg CCqgeH^ ^ "" ( 0 ^"^*^. 

and is the cubic spline interpolating operator. If q[j represents a solution 

obtained from minimizing over Q^, then there exists q eQ which minimizes J 

over Q, and a subsequence ft) •' c« } such that as -»■ <*., q^^ .»■ q . 

In the next section we present mumerical findings for double (state and 
parameter) approximation schemes such as those described here. 
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1 : Numerical Implementation and Examples. Recall from Section 2 that the 

approximating identification problem is: 

Given y, minimize J^(q) = ^ |y - 5^(t.,q)|^ over qeQ (where 

i=l 

invol>ves point evaluations, in space, of the first component of subject 
to z^(sq) satisfying the following ordinary differential equation: 

z"(t) = A'‘(q)z"(t) + p"(q)G(t;q) 

z"(0) = P"(q)zg(q). 

(We continue our discussions in terms of the transformed system (1.5) and 
criterion (1.6) even though the numerical examples summarized in this section 
involve "data" for the original system (1.1)- (1.4) used in conjunction with 
the criterion (1.7).) Since z^ e X^(q), z^ has a representation in terms of 

2N^3 

the basis elements of X^(q), z^(t;q) = I w^(t;q)B'|'(x;q). If we let [A^(q)] 

i=l 

and [f*'*] be the matrix and vector representations, respectively of A*^(q) and 
P (q)f (where f is an arbitrary function in X(q)) with respect to the basis 

elements of X^(q), and let w^(t;q) = col (w‘^^(t;q) ^ w^^(t;q) 

solves the following system of ordinary differential equations: 

w^(t;q) = [A^(q)]w^(t;q) + [G^(t;q)] 

w^(0;q) = [zQ(q)]. 

As in [5 ], this can be written more explicitly as: 
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Q^w^(t;q) = KV(t;q) + R^G(t;q) 

(5.1) 

qV(0;q) = R%(q) 

where Q and K are matrices, with elements described by (Q ). • = 

1 »J 

. = <6?. A(q)s5> , and (R^f ) . = <bJ, f> for feX(q). 
ijq ijJ 1 Jq 1 *q 

Due to the form of the B- spline basis elements we have chosen (see Section 2), 

N 

Q can be stored as a banded symmetric matrix; this banded, symmetric 

structure permits more efficient computations and requires less storage space. 

The matrix has a similar sparse (although not symmetric) structure. 

N N N 

Each element of the matrices Q and K , and of the vector R f depends 
continuously on q, therefore the representations [A^(q)3 and [f*^] are 
continuous in q. The basis elements for X (q) depend linearly on q, and hence 
are continuous in q, which implies q P*^(q)f and q T^(t;q)f (we note 
T*^(t;q) = exp(A^(q)t) since A*'*(q) is a bounded operator) are continuous mappings 
(recall this was a necessary condition in Theorem 2.1). 

We no-te that in the case where q-j and assumed to be constant, or 

to have a representation as, for example, a linear combination of spline 
elements, then the computations can be done more efficiently; in such cases, 
the numerical quadratures required to compute the inner products which form 
Q*'* and need be performed only once for each N. Then, to construct Q*'* and 
the appropriate multiples or linear combinations of these stored values are 
computed. 

Many of the computations in the software package used to generate the 
following examples were done with IMSL subroutines (for example, the optimization. 
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and the solution of the differential equation in (5.1)). Although niuch 
modification was necessary for the present application, the core of the 
package was developed by James Crowley [11]. The examples were computed 
either on an IBM VM/370, or a CDC 6600. 

The optimization is done using a Levenberg-Marquardt algorithm. For 
fixed N, each iteration in the optimization is performed as follows. Given 
q, beginning at time zero (t-j=0), a Cholesky decomposition method is used to 
solve (5.1) for w^(t;q) and w^(t^;q); this is then integrated using Gear's 
method to obtain w^(t2;q). We use the components of the vector w*^(t2;q) to 
recover z!|*(t2»q) as the linear combination of the first components of the 
basis elements. The vector C^(t2»q) is z!|'(t2;q) evaluated at each of the 
spatial observation points. Using w^(t2;q) as the initial value, (5.1) is 
solved again for te[t 2»t2], ^^{tgjq) is obtained, and this procedure is 
repeated until 5^(t^.,q) has been evaluated at all times t.; then J^'*(q) can be 
computed as the sum of the residuals, \y.- 5^(t^.,q)|^. The data {y^.} is 
read in and stored at the beginning. 

In the selection of examples to follow, the "data" has been generated 
with an independent finite difference scheme (an implicit method [17] was 
modified for our boundary conditions and the variable coefficient, q2(x)) 
applied to the model with a priori chosen "true" values q- of the parameters. 

In all examples, q-](x) is taken to be identically one (this is done to reduce 
ill-posedness, as mentioned in Section 1). We begin each example with an 
initial guess, and a value of N; we solve (ID^), to get converged values, q^ 

(these are numerical approximations (to q*^) that result from the Levenberg- 
Marquardt algorithm), which we then use as starting values for the next value of N. 
So, in Example 5.1 (below) we begin with N=4 and a guess q*^, and generate q^. We 
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-4 -8 

then start with q at N=8, and generate q . 

We remind the reader that the computations reported on below were 
carried out using "data" for the system (1.1)-(1.4) with criterion (1.7) and 
an appropriate approximate criterion for the problem. (We have also 
successfully tested the methods on similar examples with the transformed 
system (1.5) and criterion (1.6), although, of course, this is not the typical 
formulation of the inverse problem for which data will be available.) 

Example 5.1. For our first example we used "data" consisting of observations 
at x=0 and times t = .25, .5, .75, ..., 2.0. This is meant to simulate the 
situation in "surface seismic" experiments where only data at the surface are 
available. The source term was chosen as s(t;k) = qg(l-e"^^)e^^^ » a function 
which rises to a peak quickly and then gradually diminishes to zero; again 
this attempts to mimic the situation in seismic experiments. We assume van- 
ishing initial conditions and seek to estimate a constant elastic modulus qg 
as well as the boundary parameters d3» <14 and the source parameters k= (qg.qg). 

True valu'es along with our estimates are given in the results summarized in 

* 

Table 5.1, Graphs comparing the true solution at the surface u(t,0;q ) with 
the approximate solution u^(t,0;q*^) are shown in Figure 5.1. We also tested 
the method on this example using "data" for more spatial observations (data 
at x=0, .5, 1.0 and at t = .5, 1.0, 1.5) with our findings given in Table 5.2. 
Based on these computations and a number of other tests, we suggest that 
there appears to be little difficulty with our method in the case where only one 
spatial observation is available as long as a sufficient number of time 
observations are available. 


T 


34 


TABLE 5.1 


Initial 

Guess 


Converged Values 
N =4 N =8 


True 

Values 


qg = 2.0 

qg = 2.96001 

q® = 3.0001 

q2 = 3.0 

-1.0 

q3 = -1.98861 

q® = -1.99012 

qj = -2.0 

= 2.0 

qj = 0.97428 

q^ = 1.00683 

qj = 1.0 

q“= 1.5 

qg = 1.97135 

q^ = 1.99809 


q“ = -0.5 

q^ = -0.98500 

qg = -1.00506 


No. of , 

Iterations 

11 

2 


R.S.S.^ 

0.659 X 10"^ 

0.119 X 10"^ 


CPU^ 

125.363 

84.688 



^Number of iterations in the optimization algorithm. 
^Residual sum of squares = 

3 

The CPU time given in seconds. 
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TABLE 5.2 


Initial Converged Values 

Guess N = 4 N=8 Values 


q° = 2.0 


q® = 2.99378 

~ 

o 

1 

II 

O CO 

cr 

q3 = -1.92304 

qg = -2.01999 

= -2.0 

o 

II 

ro 

o 

qj = 1.01302 

q® = 1.00285 

q^ = 1.0 

4 - 

= 1.97120 

q® = 2.00578 

2.0 

q° = -0.5 


q^ = -0.99172 

q* = -1.0 

No. of 




Iterations 

12 

5 


R.S.S. 

0.235 X 10"^ 

0.441 X 10"^ 







CPU 


117.597 


147.323 
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Example 5.2. In this. example we compared the performance of our method on 
problems with "noisy data" with that on those without noise in the data. We 
used the same source term as that in Example 5.1, zero initial conditions, 
but a "true" parameterized elastic modulus E(x) = 3/2 + l/ir Arctan [q 2 i (x-q22)3* 
Data for observations at x = 0.0, 0.5, 1.0 and t = .416, .832, 1.248, 1.664, 

2.08, 2.496 were used. Results for the case of data without noise are summarized 
in Table 5.3, while findings employing data with a noise level of approximately 
3% are given in Table 5.4. In both cases, the method converges nicely but as 
one might expect, the converged values of the parameters do not agree with the 
true parameters in the case of noisy data. In Figures 5.2, 5.3, 5.4 and 5.5, 

M * 

we graphically depicted the curves for E and E in several cases. 

Example 5.3. In this example we illustrate the ideas discussed in Section 4 

regarding parameter approximation in the set of linear and cubic splines. We 

do not assume an a priori shape for the elastic modulus E(x), the "true" value 

of which is given by E (x) = 3/2 + tanh [6(x- .5)]. Rather we first search 

★ 

for E in the class of linear spline approximations to E . We then carry out 

the search using cubic splines. Initial conditions are u(0,x) = e^, 

u^(0,x) = -3e and no source term was assumed ( i.e., s = 0). Data for observations 

at 3 spatial points (x = 0.0, 0.5, 1.0) and 6 time points (t= .16, .32, ..., 1.0) 

* 

were used. Figure 5.6 depicts graphs of the true modulus E , the initial guess 
E^, and the converged estimate E^ where we used linear splines (with 4 basis 
elements -- M= 3 in the notation of Section 4) to approximate E and cubic splines 
(N=4) to approximate the state. At the same time we searched for the boundary 

•k k 

parameters qg, q^ (true values q^ = -1.0, q^ = 3.0) and obtained converged 
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= 1.0 


= 2.97352 

^21 

= 2.99994 


= 1.0 

4z 

= 0.511 15 


= 0.50053 


= -2.0 

^3 

= -0.99892 


= -1.00026 

< 

= 2.0 


= 3.05138 


= 3.01070 

^5 

= 1.0 

’5 

= 2.00322 

■4 

= 2.00056 

^6 

= -2.0 


= -1.01163 

^6 

= -1.00217 

No. 

of 





Iterations 


13 


3 

R.S. 

S. 

0.1025 X 10"^ 

0.82859 X 10"^ 

CPU 

269.696 

196.335 
















39 


Initial 

Guess 


TABLE 5.4 
(NOISY DATA) 
Converged Values 
N=4 N=8 


True 

Values 


= 1.0 q2i = 3.30536 = 3.29222 = 3.0 

q°2 = 1.0 q22 = 0.53802 q®2 = 0.53115 q*^ = 0.5 

q° = -2.0 q^ = -0.86648 q® = -0.86017 q* = -1.0 

q° = 2.0 qj = 2.99610 q® = 2.96002 q* = 3.0 

q° = 1.0 qg = 2.09207 q® = '2.09295 q* = 2.0 

0 -4 

^l6 " -2.0 qg = _i, 15602 q® =-1.15571 q* =-1.0 

No. of 

Iterations 13 2 

R.S.S. 0.6509 X lo'^ 0.476 x lO"^ 


CPU 


270.11 


136.87 
















FIGURE 5.3 
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-4 -4 

estimates qj = -1.05425, = 3.3576 with a CPU time of 38 seconds and 

_2 

R.S.S. = 0.255 X 10 . Figure 5.7 contains graphs similar to those in 

Fig. 5.6 e;<cept N=16 was used in the state approximations. Boundary 

parameter estimates corresponding to were q^^ = -1.10063, q^® = 3.07049 

with CPU time of 118 seconds and R.S.S. = 0.472 x 10"^. The error (in 

the H norm) in estimating E in each case was calculated to be |E - = .081 

and |E* - = .030. 

We carried out similar calculations for the same example in which we 
employed cubic splines (M=l in the notation of Section 4, i.e. 4 basis elements) 
for the parameter approximations. The graphs of E , E and E ° are compared 
in Figure 5.8. In this second test we did not search on the boundary 
parameters q 2 » q^ but rather held them fixed at their "true" values. The 
error at the converged parameter was |E* - E^®| = .109, with R.S.S. = 0.293 x lo"^ 
and a CPU time of 178 seconds. 



FIGURE 5.6 
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FIGURE 5.8 
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6. Concluding Remarks. 

We have presented in this paper both theoretical and numerical results 
using some of our ideas involving spline approximations for inverse or 
parameter estimation problems for hyperbolic systems. Among the novel features 
is the capability of estimating variable coefficients and boundary parameters 
with methods that are both theoretically sound and readily implementable. 

Our techniques (reported on earlier, [6]) involve the use of parameter dependent 
basis elements for the approximation subspaces in a Galerkin type semi-discrete 
scheme. 

While we have focused on 1-dimensional space domain problems here, our 
ideas. are in principle applicable to problems in 2 and 3 dimensional domains. 

We have devoted some thought to such problems in connection with use of basis 
elements that are tensor products of 1-D elements. These ideas offer some 
promise, given the parallelism that would be inherent in the resulting algo- 
rithms and given the emerging technology related to supercomputers and array 
processors. However, there are other ideas that also offer great promise; in 
particular, there are those involving spectral methods such as the tau-Legeodre for 
which we have reported preliminary findings in [4]. A fundamental -difference 
between these techniques and those proposed in this paper is that in the tau- 
Legendre one does not require the approximation subspace basis elements to 
satisfy the boundary conditions. Instead the boundary conditions are essen- 
tially imposed as side constraints adjoined to the Galerkin type differential 
equations. This can offer significant computational advantages, especially 
in higher dimensional domain problems. We are currently pursuing investigations 
of these ideas. 
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In closing we remark that the theoretical results presented above only 
guarantee convergence of subsequences to a minimizer q for J. But for 

the class of problems investigated here and for a number of other types of 
inverse problems we have studied, we have in practice only observed (numerically) 
convergence of the original sequence (;■} . This has been our experience even 
in examples with noisy data and may be due in many cases to the fact that the 
original problem of minimizing J over Q has a unique solution q . In this 
situation, elementary and quite standard arguments can be employed to actually 
establish convergence of < q > itself to q . 
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